#include <iostream>
#include <vector>
#include <list>

#include <TROOT.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1D.h>
#include <TCanvas.h>
#include <TVector3.h>

#include "GeomSvc.h"
#include "TriggerRoad.h"
#include "TriggerAnalyzer.h"

using namespace std;

int main(int argc, char* argv[])
{
  GeomSvc* p_geomSvc = GeomSvc::instance();
  p_geomSvc->init("geometry_R962");

  Double_t mass, xF, x1, x2;
  Double_t weight;
  Double_t x0, y0, z0;

  Double_t px1, py1, pz1, p1;
  Double_t px2, py2, pz2, p2;

  Int_t nHits;
  Int_t detectorIDs[1000];
  Int_t elementIDs[1000];

  TFile* dataFile = new TFile(argv[1], "READ");
  TTree* dataTree = (TTree*)dataFile->Get("save");

  dataTree->SetBranchAddress("weight", &weight);
  dataTree->SetBranchAddress("z0", &z0);
  dataTree->SetBranchAddress("mass", &mass);
  dataTree->SetBranchAddress("xF", &xF);

  dataTree->SetBranchAddress("px1", &px1);
  dataTree->SetBranchAddress("py1", &py1);
  dataTree->SetBranchAddress("pz1", &pz1);
  dataTree->SetBranchAddress("px2", &px2);
  dataTree->SetBranchAddress("py2", &py2);
  dataTree->SetBranchAddress("pz2", &pz2);
  
  dataTree->SetBranchAddress("nHits", &nHits); 
  dataTree->SetBranchAddress("detectorID", detectorIDs);
  dataTree->SetBranchAddress("elementID", elementIDs);  

  Double_t cut1 = 0.;
  Double_t cut2 = 0.;
  Double_t cut3 = 1E6;
  if(argc > 3) cut1 = atof(argv[3]);
  if(argc > 4) cut2 = atof(argv[4]);
  if(argc > 5) cut3 = atof(argv[5]);
  TriggerAnalyzer* triggerAna = new TriggerAnalyzer();
  triggerAna->init("roads_DY.root", cut1, cut2, cut3);
  triggerAna->buildTriggerTree();

  Double_t w_target = 0.;
  Double_t w_dump = 0.;
  Double_t w_target_acc = 0;
  Double_t w_dump_acc = 0.;
  Double_t w_both_acc = 0.;
  Double_t w_both_all = 0.;
  
  int nPairs;
  int p_groupIDs[100];
  int m_groupIDs[100];

  TFile* saveFile = new TFile(argv[2], "recreate");
  TTree* saveTree = dataTree->CloneTree(0);

  saveTree->Branch("nPairs", &nPairs, "nPairs/I");
  saveTree->Branch("p_groupIDs", p_groupIDs, "p_groupIDs[nPairs]/I");
  saveTree->Branch("m_groupIDs", m_groupIDs, "m_groupIDs[nPairs]/I");

  TH1D* z_all_target = new TH1D("z_all_target", "z_all_target", 25, -200., 200.);
  TH1D* z_all_dump = new TH1D("z_all_dump", "z_all_dump", 25, -200., 200.);
  TH1D* z_acc_target = new TH1D("z_acc_target", "z_acc_target", 25, -200., 200.);
  TH1D* z_acc_dump = new TH1D("z_acc_dump", "z_acc_dump", 25, -200., 200.);
  TH1D* mass_all_target = new TH1D("mass_all_target", "mass_all_target", 25, 1., 10.);
  TH1D* mass_all_dump = new TH1D("mass_all_dump", "mass_all_dump", 25, 1., 10.);
  TH1D* mass_acc_target = new TH1D("mass_acc_target", "mass_acc_target", 25, 1., 10.);
  TH1D* mass_acc_dump = new TH1D("mass_acc_dump", "mass_acc_dump", 25, 1., 10.);
  TH1D* xF_all_target = new TH1D("xF_all_target", "xF_all_target", 25, -1., 1.);
  TH1D* xF_all_dump = new TH1D("xF_all_dump", "xF_all_dump", 25, -1., 1.);
  TH1D* xF_acc_target = new TH1D("xF_acc_target", "xF_acc_target", 25, -1., 1.);
  TH1D* xF_acc_dump = new TH1D("xF_acc_dump", "xF_acc_dump", 25, -1., 1.);
  TH1D* angle_all_target = new TH1D("angle_all_target", "angle_all_target", 25, -0.05, 0.6);
  TH1D* angle_all_dump = new TH1D("angle_all_dump", "angle_all_dump", 25, -0.05, 0.6);
  TH1D* angle_acc_target = new TH1D("angle_acc_target", "angle_acc_target", 25, -0.05, 0.6);
  TH1D* angle_acc_dump = new TH1D("angle_acc_dump", "angle_acc_dump", 25, -0.05, 0.6);
  TH1D* pT_all_target = new TH1D("pT_all_target", "pT_all_target", 25, 0., 10.);
  TH1D* pT_acc_target = new TH1D("pT_acc_target", "pT_acc_target", 25, 0., 10.);
  TH1D* pT_all_dump = new TH1D("pT_all_dump", "pT_all_dump", 25, 0., 10.);
  TH1D* pT_acc_dump = new TH1D("pT_acc_dump", "pT_acc_dump", 25, 0., 10.);


  for(int i = 0; i < dataTree->GetEntries(); i++)
    {
      //cout << "=====================================" << i << endl;
      dataTree->GetEntry(i);
    
      w_both_all += weight;

      TVector3 p_vec(px1, py1, pz1);
      TVector3 m_vec(px2, py2, pz2);
      TVector3 sum = p_vec + m_vec;
      Double_t angle = p_vec.Angle(m_vec);

      if(z0 > 0) 
	{
	  z_all_dump->Fill(z0, weight);
	  mass_all_dump->Fill(mass, weight);
	  xF_all_dump->Fill(xF, weight);
	  angle_all_dump->Fill(angle, weight);
          pT_all_dump->Fill(sum.Perp(), weight);

	  if(mass > 4.) w_dump += weight;
	}
      else
	{
	  z_all_target->Fill(z0, weight);
	  mass_all_target->Fill(mass, weight);
	  xF_all_target->Fill(xF, weight);
	  angle_all_target->Fill(angle, weight);
          pT_all_target->Fill(sum.Perp(), weight);
	  
	  if(mass > 4.) w_target += weight;
	}

      if(triggerAna->acceptEvent(nHits, detectorIDs, elementIDs))
	{
	  w_both_acc += weight;

	  std::list<TriggerRoad>& p_roads_found = triggerAna->getRoadsFound(+1);
	  std::list<TriggerRoad>& m_roads_found = triggerAna->getRoadsFound(-1);
	  nPairs = 0;
	  for(std::list<TriggerRoad>::iterator iter = p_roads_found.begin(); iter != p_roads_found.end(); ++iter)
	    {
	      for(std::list<TriggerRoad>::iterator jter = m_roads_found.begin(); jter != m_roads_found.end(); ++jter)
		{
		  //do whatever you like about combination
	
		  p_groupIDs[nPairs] = iter->groupID;
		  m_groupIDs[nPairs] = jter->groupID;

		  ++nPairs;
		}
	    }
	  saveTree->Fill();
    	  
          if(z0 > 0)
	    {
    	      z_acc_dump->Fill(z0, weight);    
    	      mass_acc_dump->Fill(mass, weight);
    	      xF_acc_dump->Fill(xF, weight);
    	      angle_acc_dump->Fill(angle, weight);
              pT_acc_dump->Fill(sum.Perp(), weight);

	      if(mass > 4.) w_dump_acc += weight;
	    }	      
	  else
	    {
	      z_acc_target->Fill(z0, weight);
    	      mass_acc_target->Fill(mass, weight);
    	      xF_acc_target->Fill(xF, weight);
    	      angle_acc_target->Fill(angle, weight);
              pT_acc_target->Fill(sum.Perp(), weight);

	      if(mass > 4.) w_target_acc += weight;
	    }
	}
    }

  cout << "Total cross section before: " << w_both_all << ", after: " << w_both_acc << endl;
  cout << "Dump:   " << w_dump_acc << "  " << w_dump << "  " << w_dump_acc/w_dump << endl;
  cout << "Target: " << w_target_acc << "  " << w_target << "  " << w_target_acc/w_target << endl;
  cout << "Target/Dump before cut: " << w_target/w_dump << " ==== Target/Dump after cut: " << w_target_acc/w_dump_acc << endl;
  cout << "Target/All before cut: " << w_target/(w_dump + w_target) << " ==== Target/All after cut: " << w_target_acc/(w_dump_acc + w_target_acc) << endl;
  
  saveFile->cd();
  saveTree->Write();
  z_all_target->Write();
  mass_all_target->Write();
  xF_all_target->Write();
  angle_all_target->Write();
  z_acc_target->Write();
  mass_acc_target->Write();
  xF_acc_target->Write();
  angle_acc_target->Write();
  z_all_dump->Write();
  mass_all_dump->Write();
  xF_all_dump->Write();
  angle_all_dump->Write();
  z_acc_dump->Write();
  mass_acc_dump->Write();
  xF_acc_dump->Write();
  angle_acc_dump->Write();
  pT_all_target->Write();
  pT_acc_target->Write();
  pT_all_dump->Write();
  pT_acc_dump->Write();
  saveFile->Close();

  delete triggerAna;

  return 1;
}
